Compressive consolidation of strongly aggregated colloidal gels 
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The compressive yield stress of colloidal gels shows a highly nonlinear dependence 
on the packing fraction. We have studied continuous compression processes, and dis- 
cussed the packing fraction dependence with the particle scale rearrangements. The 
2D simulation of uniaxial compression was applied to fractal networks, and the re- 
quired compressive stresses were evaluated for a wide range of packing fractions that 
approached close packing. The compression acts to reduce the size of the characteris- 
tic structural entities (i.e. the correlation length of the structure). We observed three 
stages of compression: (I) elastic-dominant regime; (II) single-mode plastic regime, 
where the network strengths are determined by the typical length scale and the rolling 
mode; and (III) multi-mode plastic regime, where sliding mode and connection breaks 
are important. We also investigated the way of losing the fractal correlation under 
compression. It turns out that both fractal dimension Df and correlation length £ 
start to change from the early stage of compression, which is different from the usual 
assumption in theoretical models. 
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I. INTRODUCTION 



Particulate gels form in colloidal dispersions if strong short-range attractions act between 



contacting surfaces of particles. Particles once joined hardly separate (Larson, 1999 



Larson 


1999 




Witten 



and Pincus, 2004). The attractions also prevent tangential displacements of contacting parti- 



cles (Dinsmore, 2010), thus particles form sample-spanning loose networks. The remarkable 
and complex mechanical properties of a colloidal gel under stress are essentially due to the 
response of a deformable disordered network formed with all the colloidal particles. 



Under shear stress, one can sketch out deformation of a colloidal gel as follows (Nguyen 



and_BogerJ[1992J): 

• the system is elastic under very low stresses. If the stress is released, the system comes 
back to the original shape; 

• for larger stress, a colloidal gel may undergo finite deformations, though it does not 
flow. This is the plastic regime, and it gives malleability; i.e., the system mostly keeps 
the deformed shape if stress is released; 

• beyond a yield stress, the particles network is ruptured and the colloidal gel flows. 



Fluids exhibiting these behavior are called yield stress (or Bingham plastic) fluids (Nguyen 
and Boger, 1992 Denn and Bonn, 2011). Since the yield stress is the highest stress that 
can prevent flow, it can be considered as the strength of colloidal gels. Naturally, this 
strength relates to the strength of the particle connections. In a colloidal gel under shear 
stress, transmitting stresses break a portion of particle connections and the sample spanning 
network is torn apart. The yield stress of colloidal gels is higher when the packing fraction 
of particles is larger. In many cases, this variation appears to follow a power law with 



exponents between 4 and 5 (Channell and Zukoski, 1997). 



Under compressive stress, the behavior of a colloidal gel looks different, due to confine- 
ment (de Kretser, Boger, and Scales, 2003 Stickland and Buscall, 2009). Compressive 



stress appears naturally in a number of applications, for example paste drying or ceramic 
processing. Osmotic compression is a generic important case. In this kind of compression 
process, liquid removal tends to pack the particles into a smaller volume. As long as the 
applied stress is lower than the strength of the colloidal gel, it will retain the volume. If the 
stress exceeds a compressive yield stress P y , network restructuring leads to the loss of a part 



of the liquid phase, and particles reorganize into a smaller volume (Buscall et al. , 1987). The 



deformation stops when the mechanical resistance of the network equals the stress. Thus, 



unlike the shear case, deformation does not continue with a constant applied stress (Buscall 
and White, 1987). The growth of mechanical resistance, i.e. compressive consolidation, also 



shows a highly nonlinear dependence on the packing fraction, typically power-law relations 



P y ~ 4> x with exponents between 4 and 5 ( 


Miller, Melant, and Zukoski 


Boger 


1997 




Madeline et al. 


2007 




Parneix et al., 


2009 


)• 



1996 



Green and 



The highly nonlinear dependency of the compressive yield stress P y on the packing fraction 



4> relates to the heterogeneity of the network structure of particles (jBuscall et al. 
The fractal model (iBrown 



1986 




Shih et al. 


1990 





Potanin and Russel, 1996) successfully 



1987). 



captures such packing-fraction dependence of the colloidal gel's rheology by considering 
the characteristic structure of the colloidal gel. In the model, a colloidal gel is assumed 
to be an assembly of fractal floes — such structural correlation is evidenced by scattering 
observations (Dietler et al., 1986 Poon and Haw, 1997). The entire sample is filled up with 



floes of a fractal dimension Df, and the correlation length £ is equivalent or proportional to 
the average size of the floes. So, the correlation length £ can be associated with the packing 



2 



fraction by using the fractal dimension Df and the spatial dimension d: £ ~ (f)~ 1 ^ d ~ Di \ 
A higher packing fraction means that the fractal dimension D{ is the same, but a higher 
number of smaller floes fills up the space. By evaluating average mechanical responses of the 
floes, the fractal model predicts a power-law dependence on with exponent (d + D^) / (d — 
.Df), where is fractal dimension of the backbone structure. Although this theoretical 
prediction has often been cited in the literature, we should pay attention to the scope of 
the model. Since fractal correlation yields in the formation process of gels, there is no clear 
reason for Df and D hh to be kept constant during the compression process. Actually, some 
scattering observations indicate that the compression process changes the local structure 
(fractal dimension) (Madeline et al., 2007 Parneix et al., 2009). To make progress in 



macroscopic modeling of aggregated colloidal gels (Flatt and Bowen, 2006 Lester, Rudman 



and Scales, 2010), we need to settle such microstructural issues. 



There is recent simulation work to support the application of the fractal model to com- 
pression processes (Gilabert, Roux, and Castellanos, 2008). Gels which were (i) directly 



prepared to the packing fraction 0, or (ii) prepared by compression from a lower packing- 
fraction gel 0' < showed the same structural characteristics. However, the investigated 
range of packing fractions was relatively high (0.47 < < 0.78). Since the fractal model 
targets strongly flocculated gels, we should compare simulations at lower packing fractions. 
Indeed, we report here a wider range of packing fractions (0.12 < < 0.74). 

This article reports results of numerical simulation for compression of particulate networks 
to discuss the consolidation behavior with restructuring. A contact model is employed to 
simulate sticky and breakable connections between contacting particles. In granular physics, 



there are many prior studies to investigate cohesive granular system (Iwashita and Oda 



1998 Tomas, 2007). Their simulation methods are an extension of the discrete element 



method (Cundall and Strack, 1979); the rolling friction has been additionally included into 
the contact model. Though our target is a colloidal system, the spirit of the modeling is 
almost the same. In granular physics, there are also some studies on compression of cohesive 



powders ( 


Kadau et al. 


2002 




Bartels et al, 


2005 




Wolf et al. 


2005 




Gilabert, Roux, and 


Castellanos 


2007, 


2008 


), which obtained many common insights with 


our study. But this 



work is unique in its focus on strongly flocculated gels possessing long correlation lengths. 
Though we report 2D simulation results here, 2D systems are not limited as subjects of 
simulation studies, but there are some experimental investigations (Masschaele, Fransaer, 



and Vermant 2009 Cicuta and Vella] 2009 |Berhanu and Kudrolli , 2010 Park and Furst 



2011). One obvious benefit of studying 2D systems is easy recognition of the structure 



changes. Though the density-density correlation is an effective tool to capture special sorts of 
heterogeneous structures, it might not capture all features of random structures. Therefore, 
our eyes are still efficient to obtain some insights from simulation results. 



II. METHOD 

A. Initial configuration: fractal gel 

We use fractal aggregation models to generate the initial configurations. Networks of 
particles are formed in a rectangle box of sides L x and L y with periodic boundary conditions. 
We combine two processes. The first step models a formation process of fractal clusters 
in the dilute limit; the reaction-limited hierarchical cluster-cluster aggregation model was 



employed (Jullien and Botet, 1987). Each cluster consists of N Roc particles. The second step 
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is to build a network with them. N c \ of the generated clusters are randomly placed in the 
box while avoiding overlap with each other, and then they are moved with a random walk. 
Collisions of clusters always lead to coalescence, so a single cluster is formed in the periodic 
boundary space after a while (see, e.g., FIG. [7] (a), below). Thus, we can generate random 
networks of a packing fraction = TTa 2 N c \Nf[ oc /L x Ly, where a is the radius of a particle. 

One can deduce structural features of the formed networks by evaluating the density- 
density correlation function C(r) (See appendix [A]) . We prepared 10 samples of N = 5760 
particles in a box (L x ,L y ) = (300a, 500a). The packing fraction is 4> « 0.12. FIG. [I] shows 
the averages and standard deviations of their correlation functions C(r). A power-law decay 
of correlation is seen in the finite range 6 < r/a < 20. The slope indicates fractal dimension 
.Df ~ 1.56, which is almost the same fractal dimension as the clusters generated by the 
aggregation model. Though a crossover behavior is seen in the range 20 < r/a < 70, the 
correlation decays out after that. 
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FIG. 1. The density-density correlation functions were evaluated with the fractal networks of 
cf) « 0.12 prepared in the box (L x , L y ) = (300a, 500a). The averages and standard deviations were 
taken over 10 samples. The (red) dashed line indicates the slope corresponding to D{ ~ 1.56. 



B. Contact forces 



Our contact model is designed to simulate strongly flocculated colloidal gels. The contact 
model used in the original discrete element method (DEM) ( |Cundall and Stra"ck 1979) takes 
into account friction forces for the sliding mode, which is essential to simulate dense gran- 



ular systems. For cohesive granular systems, the rolling friction has been added (Iwashita 



and Oda, 1998). This extension is also required to reproduce the behavior of colloidal ag- 



gregates (Seto et ai, 2012). Our implementation is similar to the model used for cohesive 



granular systems. Since the dynamic friction, which is modeled in typical contact models, 
is unclear for colloidal systems, this part of our model is different. Though a brief summary 
will be given below, the detailed descriptions of our implementation were given in another 



2012). 



paper (Seto et al. 

We assume that interactions acting among particles are only short range forces. No 
interparticle interaction acts between two particles before getting into contact. Once the 
gap between two particles becomes zero, a cohesive force starts to act between them — our 
simulation model assumes immediate bond generations. 
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The forces and moments are assumed to Hookian: Let us consider the contact interactions 
between particles i and j. The center-to-center distance between them is r^ J ) = \r^ — r l \, 
and the gap = (n 1 '-" — 2a). The normal (separation/overlap) displacement is related 
to the normal element of the contact force, 

F§ J) = {k N 8^ + C$(-5^)(5^) 3 }n^\ (1) 

where $(x) = 1 for x > and $(x) = otherwise, n^'^ is the normal unit vector defined by 
n (M') = [ r (j) _ r w^/ r (hj) m second term of eq. (fij) is added to avoid significant overlap 
between particles. The sliding displacement d^' relates to the tangential part of the contact 
force, 

F^ 3) = k s d^\ (2) 

The sliding displacement d^' is the tangential projection of the deviation vector between 
the original contact points r^' + and r^' + a£V' i >, 

d {i ' j) = a{(£ to) - ^> j) ) - - ■ n^'W^j, (3) 

where the vectors and £W» indicate the original contact points from the respective 
centers of the particles and are fixed with them. The bending angle ip^>^ (= — sin _1 (£^'^ x 
£ )z) is related to the moment acting on the rolling mode, 

= k R aV hj) - (4) 

These elastic relations are broken when the force or moment exceed threshold values. 
We assume simple criteria given with three threshold values: the critical normal force F^ c , 
critical sliding force F$ c , and critical rolling moment Mr c . When the pulling force exceeds 
the value -Fjij • n^'^ > F^ c , the bond is broken and the two particles separate. When the 
tangential force or rolling torque exceeds the threshold valueb |i*r| > -^sc or |M R | > M Rc , 
both the sliding displacement d and rolling angle ip^> are reset to zero, which releases the 
tangential parts of the stored energy. The maximum strains are related to the threshold 
forces and moment as follows: 5 C = F^/k^, d c = F$/ks, and tp c = Mr c / '(/cro 2 ). We will test 
this model by the above-mentioned three-point bending geometry in section III A 



It must be noted that in our simple model the bond strength of the tangential mode does 
not depend on the normal load. However, the enlarged contact area due to the normal load 



may change the tangential strengths (Dominik and Tielens, 1995). Our neglect of this effect 



can affect simulations at higher packing fractions. 



C. Equation of motion 

A particle i may be in contact with other particles j. The total contact force and torque 
acting on the particle i can be written as follows: 



j , \ (5) 

jf = {Mi' 3) + on^ x F^ 
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Particles are immersed in a viscous fluid, so that hydrodynamic interactions also act 
on them. For micro or nano-sized particles, the particle Reynolds number is extremely 
small, so the Stokes equations should be used to evaluate the hydrodynamic interactions 



among many particles (Kim and Karrila, 1991). However, we do not go into the detailed 



hydrodynamics here. Hydrodynamic interactions in the Stokes regime are scaled with the 
velocities of particles. If the slow compression limit is considered, flows induced by particle 
movements are negligible in comparison with the contact forces. They just work as energy 
dissipators. As an approximation, the Stokes formula is used in this work; the force and 
torque acting on a particle i are proportional to the translational velocity and angular 
velocity respectively: 



rp{i) _ 



-d7TT] av 



in 



(6) 



where 7/0 is the viscosity of the solvent. 

Since the interparticle interactions are strong enough in comparison with thermal agita- 
tions, Brownian forces are also neglected in this work. Consequently, Newton's equations of 
motions are given as the sum of the contact forces ^ and hydrodynamic interactions (|6]), 



dt c ' 



2 
-? 

5 



(7) 



Since particle inertia is very small for colloidal systems, it is a good approximation to 



neglect the inertia terms (the left hand sides) and solve the force balance equations (Ball 



and Melrose, 1997). As explained later (sect. HE), however, we took a different strategy to 



find the static equilibriums. 



D. Uni-axial compression and static equilibrium 

The compression is externally controlled. For simplicity, we focus on examining the uni- 
axial compression along the y axis in this work. The periodic boundaries at y — and 



y = L y of the simulation box introduced in section II A are replaced by two flat walls. These 
walls can move along the y-axis like pistons, and they are assumed to be semipermeable; 
liquid passes through them, while particles do not. When particles get in contact with the 
walls, they are rigidly fixed there; particles contacting the wall move together as one object. 
The periodic boundaries along the x-axis are still used during the compression simulation. 

The stress-controlled compression, in which an external stress is applied to cause the 
compression, is similar to typical experimental situations, such as osmotic compression. 
However, due to the higher efficiency of the simulation, we use strain-controlled compression. 
This also allows us to have better statistics in our data. This work investigates the slow 
compression limit; the time scale of compression is much longer than the relaxation time of 
particle motions. In this case, the results obtained with the strain-controlled compression 
are essentially equivalent to the ones with the stress-controlled compression. 

The compressive stress balancing the network resistance is evaluated from the total forces 
acting on the walls from particles. Since particles directly connecting to the walls (called 
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wall particles) no longer have any mechanical freedom, the forces between wall particles i 
and particles j contacting them are summed up as follows: 

i j 

When we refer to "compressive stress" in this 2D simulation work, our simulation box is 
assumed to have a thickness of particle diameter 2a. Therefore, the compressive stresses are 
given as P top = F waa /(2aL x ) and P hot = -F wa31 /(2aL x ). 

Our compression simulations are performed in a stepwise manner in order to obtain 
numbers of equilibriums by a single simulation run. Compression steps are given by two 
increment ratios wq and w; the initial compression from 0o to <pi is given as <fi\ = wo4>o, and 
the k-ih packing fraction as <pk = w k ~ 1 <pi. The increment ratio w is determined with the 
number of steps k' from (pi to max : w = ((p max / '(pi) l ^ k ~ l \ The strain-controlled compression 
is simply realized by translations of the walls with the velocity f wa ii, so that the packing 
fraction increases continuously. 

When it reaches a target packing fraction the translations of the walls are suspended 
to await the equilibrium. Usually, the excess stress due to the quickness of the compression 
is relaxed by spreading the deformations into the entire system. The following criteria are 
examined every m exam time steps of the simulation: 

• The relative difference between the top and bottom stresses is smaller than a threshold 
value A: |P top - P ho t\/P < A where P is the average, P = (P top + Pbot)/2. 

• The relative change of the compressive stress for the last m steps is smaller than a 
threshold value B: \P — P'\/P < B, where P' is the one of m exam steps before. 

• No bonds were broken for the last m exam steps. 

Thus, the compressive stress P& balancing the network of packing fraction (pk can be deter- 
mined. Once Pfc is determined, the compression continues again for the next target (pk+i- 



E. Non real-time simulation for quasi-static simulation 

In order to reduce finite size effects, the simulation box needs to involve a sufficient 
number of typical structures. The compressive stress to balance sparse networks is very 
low, and the most of particles are slightly stressed. In this case, the relaxation time is long. 
Usually, the overdamped equations of motion, which are given by dropping the right hand 
sides of eqs. (|7|) due to the negligible inertia, are used to simulate colloidal suspensions. 
However, we have taken a different approach for efficiency. As restricting our scope in the 
slow compression, we may be allowed to optimize some parameters for a faster simulation. 
In the slow compression limit, velocities of particles are also slow. In this case, the contact 
forces dominate the displacements of particles. This is why a reduced viscosity rf Q still leads 
to the same equilibrium configurations with a shorter time. The following viscosity gives 
the critical damped motion in terms of the rolling mode. 



Ana 

Though it is not simple problem to find the best value, we used r]' = 0.1?7o Cd ^ in the equation 
of motion ([7]). This arbitrary selection of the viscosity is justified only if the compression is 



considered slow enough. This point will be checked in sect. HI C 
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III. RESULTS AND DISCUSSION 



A. Parameters for contact model 



The introduced 2D contact model (see sect. II B) is relatively simple, but it still has 7 
parameters to be given. Traditional experimental measurements to characterize the particle- 
particle adhesion determine the pull-off force (Larsen, 1958), which corresponds to F^ c in 
this work. In recent years the atomic force microscope became a reliable method, which 



provides more precise data for normal force-displacement relationships (Hodges et al 
However 



2002). 



there are fewer experimental measurements to determine the parameters of the 

The experimental measurement of 



other modes (Heim et al. , 1999 Ecke et al, 2001). 



the three-point bending test by using optical tweezers is remarkable (Pantina and Furst 



2005 Furst and Pantina, 2007), which provides parameters for the rolling mode, Mr c and 



k-R (or (p c ). In any case, all parameters for a single system have never been determined so 
far. Thus our knowledge of the contact forces is still limited, so that we need to set the 
parameters by hand. But, we may have a certain level of confidence by reproducing a similar 
bending behavior with the three-point bending measurements. Our selected parameters are 
intended to imitate the observed mechanical behavior: a small extent of elastic deformation 
and irreversible reorganization due to a critical moment. 

We have selected two sets of the parameters called bond 1 and bond 2, which are shown 
in TABLE [Tj Bond 1 has the same strengths for the three modes, while bond 2 is 5 times 
stronger for the normal and sliding mode than the rolling mode. Both have the same 
maximum strains; elastic deformations give 5% of the radius. 



TABLE I. Parameters for the contact model. C is the parameter in eq. ([I]) to avoid significant 
overlap by compression. 









bond 1 


bond 2 


Normal strength 




F Nc /F 


1 


5 


Sliding strength 




F Sc /F 


1 


5 


Rolling strength 




M Rc /F 


1 


1 


Normal strain 




5 C / a 


0.05 


0.05 


Sliding strain 




d c /a 


0.05 


0.05 


Rolling strain 




fc 


0.05 


0.05 


Parameter in eq. 


(il 


C(a 3 /F ) 


8 x 10 4 


8 x 10 4 



The simulation for the three-point bending test consists of the following three steps: 

i. Prepare a linear chain under no stress. 

ii. Apply the forces (F ex , 0, 0) at the central particle and (— F ex /2, 0, 0) at the both ends, 
and wait for equilibrium. 

iii. Stop applying the forces, and wait for equilibrium. 

The three snapshots i hi in FIG. [2] (a) and (b) show the configurations of the respective 
equilibriums of the bond 1 simulation. Though the data is not shown here, the results of 
bond 2 are similar. The two results were obtained with two slightly different forces: (a) 
Fex/Fo — 0.22, (b) F ex / F = 0.23. If the external force did not reach the threshold, the 
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initial linear shape remained after the test (see (a) iii). On the other hand, if the bond stress 
exceeded the threshold value, some irreversible rearrangement took place. As a result, the 
initial shape was not recovered (see (b)iii). In the latter case, the bending deformation 
proceeds until the moment acting on the broken bond becomes smaller than the threshold 
due to the increase of the bending angle. This is why the small increment of the applied 
forces leads to the visible difference of the deformation (see ii of (a) and (b) in FIG. [2]). The 
slight bending deformation seen by (a) ii is close to the limit of elastic deformation of this 
cluster. 

Thus, the parameters were chosen to simulate some brittle behavior of colloidal aggre- 
gates; the local bushy structures remain unless bond breakups take place. In this test, the 
rolling mode of the contact model is responsible for the cluster's deformation. 




FIG. 2. The three-point bending test for a linear aggregate consists of the three steps i-iii. When 
a smaller external force (F ex = 0.22Fo) is applied, the elastic deformation is seen (a), and when it 
is slightly larger (F ex = 0.23Fq), the bond breakups cause the plastic deformation (b). 



B. Parameters of the compression simulation 

The initial configurations are fractal-like with a packing fraction 0o = Nit/L x L y w 0.12, 
where (L x ,L y ) = (300a, 500a) and N = 5760 for bond 1 simulation. Since the calculation 
cost of the bond 2 simulation is higher, smaller system size, (L x ,L y ) = (200a, 500a) and 
N = 3840, is simulated. The first target packing fraction is set to <pi = (4/3)0o ~ 0.16. The 
increment ratio and the number of compression steps are w = 1.1 and k' = 17, respectively. 
The compression reaches to </> max ~ 0.74. The relaxation time depends on the parameters 
of the bonds. To account for this, the wall velocities are set to v wau n/vo = 5 x 10~ 4 and 
Vvtaxx/vo = 1 x 10~ 4 for the bond 1 and bond 2 simulations, respectively, where Vq = a^fk^Jm. 
In the equilibrium checks after every m exam = 10 4 time integration steps, A = 10 -2 is used 
for the both simulations, and B = 10~ 3 and 10~ 4 for bond 1 and for bond 2 simulations, 
respectively. 



C. Uniformity of compression 



As explained in sect. HE it is useful to introduce artificial parameters to speed up 



the compression simulation. Therefore we started to check that our method reproduced the 
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quasi-static compression as expected. If compression is too fast, the regions near the wall are 
easily crushed. So, we examined the uniformity of the deformations along the compression 
axis. The uniformity of compression can be seen via the y dependence of the packing fraction 
(ft(y), which is defined as the packing fraction in the sliced ranges: y — 5L/2 < y' < y + SL/2 
(where 5L = L y /20 was chosen). Due to the heterogeneous structure and finite system size, 
the uniformity cannot be judged with a single compression simulation. This is why we took 
the averages over 10 simulation runs with different initial configurations. 

FIG. [3] shows the results for the bond 1 simulations, where every second compression step 
is plotted. Although we can see an undesirable drop at the middle of the first compression 
((ft = 0.16), the compression is satisfactorily uniform along y-axis after that. The steep drops 
at the both edges relate to the integrity of the walls; particles connecting to the walls are 
not allowed to move (see section II D). Equivalent results were also obtained with the bond 



2 simulation with a slower compaction speed (see sect. Ill B ) 




-200 -100 

y/a 



100 200 



FIG. 3. The y dependence of packing fractions <ft(y) are shown for every second compression step 
((pi, i = 1, 3, 5 ... ). The dotted line indicates the initial configuration. 



D. Data of compressive consolidation 

Once the system reaches a static equilibrium under a compressive stress P, further com- 
pression requires a finite increment: P — > P + AP. Therefore, the balancing stress P can 



be considered as the strength of the gel network, i.e., the compressive yield stress P y (Bus- 



call and White, 1987). A simulation run of the stepwise compression (sect. HE) finds a 
series of P y vs. (ft relations. FIG. [4] plots the results of the bond 1 and bond 2 simulations, 
where the unit of the stress is Pq = Fo/a 2 . The compression from (ft « 0.16 to (ft ~ 0.74 
shows wide ranges of P y (about four decades). These results indicate significantly rapid 
compressive-strain hardening. 

Simulation provides various data, which allow us to understand the compression. Some 
of data indicate that there are three regimes in the compression process. Let us see them 
below. 

Bond creation Bond creation is essential to understand the irreversible consolidation of 
the compression. Due to heterogeneous structure of the networks, the compressive deforma- 
tions are not affine. Consequently, two separated particles may come in contact and newly 
bond with each other. The mean contact number per particle reflects new bond creations. 
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FIG. 4. The compressive yield stress Py(4>) is evaluated by the simulations of stepwise compression. 
Its averages and standard deviations taken over 10 runs are shown. The red squares represent the 
result of the bond 1 simulations, and the blue circles of the bond 2 simulations. The dotted curve 
shows the fitting function of eq. (11), and (solid and dashed) straight lines show the power-law 
functions (12). The arrows indicate ranges of the three regimes: (I) elastic-dominant regime, (II) 
single-mode plastic regime, and (III) multi-mode plastic regime. The average exponents A are 
written on the figure. 
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FIG. 5. (a) The mean contact number per particles increases as the compression proceeds. The 
standard deviations are taken over the averages of 10 simulations, (b) The successive increments 
of mean bond stored and mean dissipated energy rates are shown, where Eq = k-Qa 2 (pl/2. The 
rates are obtained by normalizing with the compressive strain A log <fi. The averages and standard 
deviations are taken over 10 runs of the bond 1 simulation. 



The result of bond 1 simulation is shown in FIG. [5] (a). The compression continuously 
changes the value from 2 to about 4 with increasing rate, so we do not see distinct regimes 
in this. 

Stored and dissipated energies The bonds have finite strength, and they can be broken 
up or irreversibly reorganized by certain applied stress. Such events cause the plastic defor- 
mation of the networks. The translations of the wall by external stress inject energy into the 
system. Part of the energy are stored in contact bonds first and dissipated as bond ruptures. 
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We can see this energy flow in the successive increments of the total stored and dissipated 
energies for each compression step (FIG. [5] (b)). The total dissipated energy exceeds the 
total stored energy at the early stage, ~ 0.2. Accordingly, we can identify the elastic 
dominant regime below and plastic regimes above. 

Bond rupture The accumulated numbers of the bond ruptures Nr™ were recorded as 
distinguishing the stress modes causing the breakup, where the mode indicates 'connection 
break', 'sliding', or 'rolling'. The rupture rates Trup de ^ can be defined as the number of 
breakups normalized with the total particle number and relative compressive strain: 



i A A T (mode) 
y{mode) _ 1 aiv rup 

mp ~ N dln</> 



(10) 



FIG. [6] shows the results of the both bond 1 and bond 2 simulations. At lower packing 
fractions, only rupture events occur due to rolling. This indicates that compression of loose 
networks is achieved by bending deformations of particulate chains. When the compaction 
reaches a certain point, further compression requires sliding or separation of contacting 
particles. The onset of this mixed-rupture mode depends on the parameters of the contact 
forces. The connection breaks start between <ft = 0.31 and 0.34 in bond 1 simulation, and 
between cf) = 0.42 and 0.46 in bond 2 simulation. They split the plastic deformation into 
two parts: single-mode plastic regime below and multi-mode plastic regime above. 



(a) Bond 1 — (P)-+- ^ — (in) 

1 r 



100 
10 



c 1 



0.1 
0.01 



i 1 1 r 

x — i connection break 



sliding 
rolling 



m 


» * i 


© 

m 


i 




I 


I 




I 


S 


, f 


s 

1 


0.3 


0.4 








S * © « 9 



(b) Bond 2 

100 



10 

a 

2 1 

4 

0.1 
0.01 



-(H)- 



-(m)_ 



— i 1 1 1 

-x — i connection break 



sliding 
rolling 



i , 

i 



0.1 0.15 0.2 0.3 0.4 0.6 0.8 1 



FIG. 6. Bond rupture rates r rup are compared for two types of bond: (a) bond 1 and (b) bond 2. 
The bond ruptures are recorded by the causing stress. The rates of the separation due to normal 
forces are plotted by squares, the ones of the regeneration due to sliding forces or rolling moments 
by triangles and circles, respectively. 



E. Three stages of compression 

The date shown in the previous section indicate that three distinct behaviors can be 
recognized in the compression process. Let us discuss them one by one here. 

(I) Elastic dominant regime There is no percolation path in the 10 initial configurations 
of 0o ~ 0.12 (see FIG. [7] (a)). Hence, no compressive stress is required to compress them 
in the slow compression limit. After some compaction, the first spanning path appears, and 
the network begins to resist the compressive stress. At the gel point, a single-linking chain 
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(percolation path) resists the external stress between top and bottom. This system-spanning 
cluster is deformed as a whole, so that strains of single bonds can be very small. Therefore, 
the dominant response to the external stress is elastic deformation; the work by compression 
is stored in the contact bonds (see FIG. [5] (b)). 

In this range, the power-law form expanding at the gel point g is suitable to capture the 



P y -(j) curve (de Gennes, 1976): 



Py/P = C\<t> 



for 



> 



(11) 



This function allows us to estimate the gel point <j) g by extrapolation. In FIG. |4j the dotted 
curve shows the fitting result with the function (11) for the bond 1 simulation data of 
4> g < (p < 0.2. The gel point was found as (fi g ~ 0.15. The evaluated exponent A' ~ 1.5(^ 1) 
represents some nonlinearity. The compression induces some new contacts, so that the 
numbers of percolation paths are increased. The rise of resistance due to such increasing 
paths can be a plausible interpretation of the nonlinearity in this elastic-dominant regime. 

It must be noted that our simulation results in this range might have large statistical 
fluctuations. The typical length scale of the structure is comparable to the system size 
(compare the drawn circle and the box size in FIG. [7] (a)), and consequently only a few 
paths enter the averaging. This is seen in the error bars in FIG. |4j 

(II) Single-mode plastic regime By considering two data (FIG. [5] (b) and FIG. [6]), we 
identified the single-mode plastic regime in the ranges: 0.2 < <p < 0.31 for bond 1 simulation 
and 0.2 < <fi < 0.42 for bond 2 simulation. 

This ^-dependence of the consolidation behavior has been discussed with a scaling con- 
cept (Buscall et al, 1988 Shih et al. 1990). The power-law relationships of P y and <fi have 



et al. 



been also observed in experiments (Buscall et al., 1987 Madeline et al., 2007 Parneix 



2009J): 

P y (<P)/P = C<P X . (12) 
A similar power-law behavior can be seen in our simulation results (FIG. [4]). The least 
squares fitting gives the almost same exponents: ~ 5.8 ± 0.5 and A^ 
bond 1 and bond 2 results, respectively. 



5.8 ± 0.4, for 



Let us look more closely at this regime. FIG. [7] (b) shows the configuration of particles 
at the beginning of this range (<p ~ 0.21) in a bond 1 simulation. At the gel point, a single 
linking path spanned the entire box. In contrast, multi-linking structures span the system 
here. Since the maximum bending moment of a random-walk chain is proportional to the 
inverse of the chain size, the supportable loads of larger loops are smaller. The variety of 
the mesh sizes indicates a heterogeneity of the local strengths. Though large empty spaces 
remain, the relatively dense meshes fill the space from the bottom to the top (FIG. [7] (b)). 
Further compression requires destroying a part of them. As shown in FIG. |6j the rolling- 
mode ruptures are the only rupture events seen in this range. Therefore, the size of the 
responsible mesh structures determines the compressive yield stress P y . 

Thus, in the single-mode plastic regime, the length scale is essential; the compression 
breaks larger and weaker structures into smaller and more robust structures. Fractal model 



was employed to give the length scale from the packing fraction (Brown, 1986 Shih et al. 
1990) (see sect. Ill F). According to the nature of a fractal, a portion of a fractal cluster is also 
a small fractal cluster. If the compression of fractal networks causes only fragmentation of 
fractal clusters into smaller pieces, it can be expected that the correlation length is decreased 
while keeping the same fractal dimension. This optimistic expectation leads us to employ 
the fractal model for the compression problem. We will check this point in the next section. 
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FIG. 7. Snapshots of the equilibrium configurations under uniaxial compression, 
configuration (a) is a prepared sample in sect. II A, whose packing fraction is (ft f 



The initial 
» 0.12. The 
b ps 0.21) and 
The radius of 



configurations (b) and (c) show the beginning of the single-mode plastic regime 
the multi-mode plastic regime ((ft ~ 0.34), respectively, and (d) shows (ft ~ 0-61 
red circles in (a)-(c) shows the correlation length £, which will be evaluated in the later section 
(sect.|mFft. 



(Ill) Multi-mode plastic regime The onset of connection breaks indicates that the com- 
pression enters a new regime. This transition point depends on the parameters of the bonds. 
We consider the data above (ft ~ 0.34 and (ft ~ 0.46 in bond 1 and 2 simulations as the multi- 
mode plastic regime. In these ranges, the compression curves look to follow a power law. 
However, the evaluated exponents are totally different with each other: Ajjj ~ 4.3 ± 0.2 and 
A§S « 5.8 ± 0.4. 

FIG. [71(c) depicts the beginning of this regime of one bond 1 simulation ((ft 0.34). There 
are still large loops remaining, but they look like voids due to the denser surroundings. Here, 
we can easily find dense meshes and trimers (a trimer can be considered as the building block 
of a close packing). If they are next to each other, the entire object can be considered as 
a lump. The same scenario of the single-mode plastic regime seems to work; denser mesh 
domains are piled up over the entire system. However, due to the finite size of particles, the 
mechanical response of denser meshes is not the same as looser ones. The crush of them 
involves sliding ruptures and connection breaks, which are seen in FIG. |6j Since the criteria 
of destruction includes the all three modes, the compression curve depends on the details of 
the contact force. 

This regime terminates due to the close packing, where rigidity of particles prevents 
reorganizations of particles. Since our simplified contact model is not accurate for that, we 
cannot specify the upper end of the range. 
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F. Fractal correlation under compression 



During compression, separated particles can newly come into contact. Consequently, 



mesh size decreases as is increased by compression. If we employ the fractal model (Brown 



1986 Shih et al. , 1990) to explain the 0-dependence of P y , the structure change needs to 
follow a special manner: the fractal dimension Df characterizing the local structures does 
not depend on the packing fraction 0, and the correlation length £ is a decreasing function 
of the packing fraction 0, £ ~ 0-i/( rf - £) f). 

The fractal dimension D{ and correlation length £ can be estimated from the density- 
density correlation function C(r) (see Appendix[A]). FIG.[8]shows averages of the correlation 
functions C(r) of every second compression step, fa (i — 1, 3, 5, . . . ), of the bond 1 simulation 
results. The results of the bond 2 simulations were almost the same. The correlation of the 
density is mostly decreasing functions of the two-point distance r. The correlation of the 
configuration before compression is the largest and reaches the farthest distance (dashed 
line). The compression reduces this correlation in a continuous manner. 
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FIG. 8. The density-density correlation functions C(r) of compressed networks are shown. The 
dashed line shows the initial configuration 0o, and the colored solid lines the every second com- 
pression step 4>i (i = 1,3,5,... ). They are averaged values over 10 runs of simulations. 



The fractal correlation appears as a straight slope in the log-log plot of C(r), i.e., 
C{r) ~ r -( d - D i). Though the straight lines can be recognized in the initial configuration, 
the compressed networks show slightly curved lines. In order to make them easier to see, we 
introduce the fractal dimension profile -Df(r) from the local slope of the correlation function 
at each distance r: 

dlnC(r) 



DAr) =d + 



dlnr 



(13) 



As mentioned in sect. II A| the initial configurations 

1.56, 



The results are shown in FIG. |9j (a 

(dashed line) consist of three parts: a fractal plateau (6 < r/a < 20) indicating D* { 
crossover range (20 < r/a < 70), and 2D-like uniform network (r/a > 70). Thus, though the 
recognizable fractal plateau is limited, we can confirm that the -Df(r) profiles of the initial 
configurations agree with the expected behavior as fractal networks. 

Once the networks are compressed, the fractal correlations are affected immediately. 
Indeed, the first compression step of our simulation (0 w 0.16) shows the entire change 
of .Df(r): a narrower plateau indicating a rise of the fractal dimension, and shifts of the 
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FIG. 9. (a) The fractal dimension profiles Df(r) evaluated by eq. (13) are plotted. The dashed 
line shows the initial configuration (fto , and the colored solid lines the every second compression step 
(pi (i = 1, 3, 5, . . . ). The numerical values indicate the packing fractions, (b) The 0-dependence of 
the correlation length £. 



crossover and uniform ranges to shorter lengths. These results suggest that not only long- 
distance correlations, but short-distance correlations can be affected by the compressive 
deformation at the early stage. The same tendency is seen in the further compression. 

Though the profiles -Df(r) for the compressed networks indicate that mono-fractal as- 
sumption does not work, it may be worth determining the <fi dependence of the correlation 
length. Here, we consider here the extent of the correlation, i.e., the length scale which 
identifies the network as 2D-like uniform network, as correlation length £. Since the profiles 
.Df (r) are not so smooth due to the finite system size of the simulation, we need to introduce 
a threshold dimension d! < d to judge the uniformity. So, an estimation of the correlation 
length £ is given as the smallest length r satisfying Df(r) = d' . Since the correlation func- 
tions of close distances oscillate due to the contacting structures of spherical particles, we 
cannot determine the fractal correlation for closer distances r/a < 6. Though this cutoff d' 
causes some underestimation of the correlation length, the obtained values are shown here. 
The results with d! = 1.96 are shown in FIG. [9] (b). According to this determination, the 
correlation length shrinks from £ w 62a at <p « 0.12 to £ ~ 8.3a at <fi w 0.41. This 
dependence follows a power-law with an exponent —2, which is a little bit slower than the 
fractal assumption: — l/(d — Df) ~ —2.2. The evaluated correlation lengths are also shown 
as the radii of the circles drawn in the configurations of particles (FIG. [7] (a)-(c)). 



IV. CONCLUSION 



We focused on investigating compression of colloidal gels. Inspired by scaling theory 
for polymer gels, the fractal model has been introduced to explain the packing fraction 
dependence of colloidal gels' rheology (Brown, 1986 Shih et al, 1990), and referred for 



the compression process (Buscall et al., 1988). The model assumes that the structure of a 



colloidal gel is fractal network. However, the simple question whether a colloidal gel under 
compression is fractal network or not, has left unanswered. We have aimed to answer this 
question by introducing a particle-scale simulation model in the limit of strongly flocculated 
gels. In this case, the cohesive force among contacting particles is assumed to be so strong 
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that no free rolling is allowed and Brownian motion is negligible. 

Our simulation reproduced the reasonable consolidation behaviors as vs. P y relation- 
ships in the case of quasi-static compression. By monitoring various particle-scale informa- 
tion, i.e., mean contact number, stored/dissipated energy, and number of bond ruptures, we 
identified the three distinct regimes: (I) elastic dominant regime, (II) single-mode plastic 
regime, and (III) multi-mode plastic regime. At the gel point, there are many "dead end" 
branches, where the particles in such branches do not contribute to the strength of the 
network P y . Compression makes separated particles come into contact, so that the mesh 
size becomes smaller and smaller. This explanation is already known and qualitatively true 
for the most part. Indeed, such shrinkage of the typical length scale causes a transition of 
the local mechanical response from (II) to (III) in the plastic regimes. While the length 
scale is large enough, the rolling mode is the responsible for the deformation. At a certain 
length scale, the sliding and normal modes become more important for further compression. 
Thus, we confirmed that a single regime could never cover the wide range of packing frac- 
tions. The simulation results also provided a insight that heterogeneity in the local mesh 
sizes plays some roles in the way of reorganization under compression and overall network 
strength. Mechanical correlations (connectivity) determine which parts of the network are 
destroyed under an applied stress. If a percolation path consisting of robust domains is 
built, they support a large fraction of the applied stress. In this case, the weakest domain in 
the path is most probable to be destroyed. This restructuring yields a smaller and robuster 
domain, which results in the reinforcement of the percolation path. However, such mecha- 
nism is subtle because compression is not continuing deformation and tends to reduce the 
heterogeneity. 

In order to see whether the fractal assumption remains or not in colloidal gels under 
compression, we evaluated the density-density correlation C(r) for the simulation results. 
The evaluated fractal dimension profile -Df(r) says that compression affects the both local 
structure and correlation length in the early stage. Thus, our results are contradictory to the 
application of the fractal model to compression processes. However, the fractal dimension 
profiles Df(r) show a systematic change in the single-mode plastic regime. This result seems 
to indicate required ingredients to upgrade the fractal model for compression processes. 
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Appendix A: Density-density correlation function 

In order to study colloidal gels, we need to have some tools to characterize random 
structured systems. The density-density correlation function C(dr) is used for this purpose. 
The local density function p(r) is defined as follows: 




1, r is in the solid phase, 
0, r is in the liquid phase. 



(Al) 
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By taking the average over the entire space of the system, the density-density correlation 
function C(dr) is defined as follows: 

(P[r)) 

When the system is isotropic, the correlation is given in terms of the distance, r = \dr\, 
which can be obtained by taking averages over all directions: 

C(r) = (C(dr)) r=ldrl (A3) 

The correlation function C(r) can be regarded as a conditional probability to find particles. 
Since it is normalized by the entire probability, the value 1 means no correlation. 

It is known that strongly flocculated gels have structures involving fractal correlations. 
The local structures are equivalent to fractal floes. For fractal floes, the relation between 
the number of particles N and the radius of gyration R g follows a power-law: 

iV ~ Rg ! . (A4) 

The exponent Df is called fractal dimension. Therefore, the average packing-fraction profile 
of the clusters is given as follows: 



r -^ Dl \ (A5) 
where r is the distance from the center of mass. Since fractal gels can be considered as 



cluster-filling networks, the local structure should coincide with this density profile (A5). 
Consequently, the density-density correlation function C(r) has the same functional form: 
C{f) ~ <f>{r). In a fractal gel, the correlation reaches a finite distance, which is called the 
correlation length £. For r > £, the correlation does not remain: 




C(r)~l ' (A6) 

I 1, for r > 4. 

If the local densities are evaluated in boxes of size £, the system should look uniform. 
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